{
"cells": [
{
"cell_type": "markdown",
"id": "analyzed-founder",
"metadata": {},
"source": [
"This tutorial presents another example of solving ordinary differential equations using *odeint()*."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "unusual-income",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.integrate import odeint"
]
},
{
"cell_type": "markdown",
"id": "handmade-corporation",
"metadata": {},
"source": [
"The differential equation that we will attempt to solve is:\n",
"\n",
"$\\ddot{z}=-\\dfrac{1}{z}\\left(\\dot{z}^2+b\\dot{z}+gz-gh\\right)$.\n",
"\n",
"In this expression, $z$ the height of water inside a straw with one end partially submerged in a cup of water. Oscillations of the water level are induced by an initial pressure difference between the inside of the straw and the surface of the water in the cup (at one atmosphere). For more details, see R. P. Smith and E. H. Matlis, *American Journal of Physics* **87**, 433 (2019).\n",
"\n",
"> $g = 9.81~\\mathrm{m}/\\mathrm{s}^2$ \n",
" $h$ is the depth that the straw is submerged in the water \n",
" $b$ is a drag coefficient"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "frank-investor",
"metadata": {},
"outputs": [],
"source": [
"g = 9.81 # m/s^2\n",
"h = 0.10 # m\n",
"b = 0.25 # m/s"
]
},
{
"cell_type": "markdown",
"id": "documented-latvia",
"metadata": {},
"source": [
"The strategy to solve a second-order differential equation using *odeint()* is to write the equation as a system of two first-order equations. This is achieved by first writing $x[1] = \\dot{z}$ and $x[0] = z$. In that case, or original second-order equation can be expressed as:\n",
"\n",
"$\\ddot{z}=\\dot{x}[1]=-\\dfrac{1}{x[0]}\\left(x[1]^2+bx[1]+gx[0]-gh\\right)$.\n",
"\n",
"One of our first-order equations is the expression above and the other is simply $\\dot{z}=x[1]$. \n",
"\n",
"- Now, we define a function that returns $\\dot{z}$ and $\\ddot{z}$ (in that order)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "modular-payroll",
"metadata": {},
"outputs": [],
"source": [
"def z_derivatives(x, t):\n",
" return [x[1], -(1/x[0])*(x[1]**2 + b*x[1] + g*x[0] - g*h)]"
]
},
{
"cell_type": "markdown",
"id": "patent-costs",
"metadata": {},
"source": [
"- Next, we create an array of the desired times for the solution."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "ancient-edinburgh",
"metadata": {},
"outputs": [],
"source": [
"time = np.arange(0, 3, 1e-3)"
]
},
{
"cell_type": "markdown",
"id": "tough-landing",
"metadata": {},
"source": [
"- Here is the call to *odeint()*. This time we need to pass an array of initial conditions, the first is for $z(0)$ and the second is for $\\dot{z}(0)$. The .T is necessary so that we can separately unpack the solutions as $z(t)$ and $\\dot{z}(t)$ (not intuitive, in my opinion)."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "cardiac-baseline",
"metadata": {},
"outputs": [],
"source": [
"position, velocity = odeint(z_derivatives, [2e-3, 0], time).T "
]
},
{
"cell_type": "markdown",
"id": "intellectual-approach",
"metadata": {},
"source": [
"- A plot of the fluid height (position) inside the straw as a function of time."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "junior-driving",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABUKElEQVR4nO3dd3yV5f3/8dedk5OEQBIIOxAIshFEBFHEhUjco9aJg2pbtVoXrbtWsNbR/hx14fyq1SK2VawVB7GylCEbFUWmICvMJJB1kty/Py5OTiIBzknuc+4z3s/HIw+uc+c+9/3h4ubkk2tatm3biIiIiMSoJLcDEBEREWkKJTMiIiIS05TMiIiISExTMiMiIiIxTcmMiIiIxDQlMyIiIhLTlMyIiIhITEt2O4Bwq6mpYdOmTWRkZGBZltvhiIiISBBs26akpIScnBySkg7e9hL3ycymTZvIzc11OwwRERFphA0bNtC5c+eDnhP3yUxGRgZgKiMzM9PRa/t8PqZOnUp+fj5er9fRa8cb1VXwVFfBU10FT3UVPNVVaMJVX8XFxeTm5tb+HD+YuE9m/F1LmZmZYUlm0tPTyczM1AN/CKqr4Kmugqe6Cp7qKniqq9CEu76CGSKiAcAiIiIS05TMiIiISExTMiMiIiIxTcmMiIiIxDQlMyIiIhLTlMyIiIhITFMyIyIiIjFNyYyIiIjENCUzIiIiEtOUzIiIiEhMUzIjIiIiMU3JjIiIiMS0uN9oUn7Ctwe2fAp710FKNrQfAc1z3Y5KRESk0ZTMJArbhu+fga/uh8pdgeNWEnS7Co56AlJauhaeiIhIY6mbKRHU+GDOVbDw5vqJDIBdA2teg4+Ogj1rXQlPRESkKZTMxDvbhvk3wro3A8e6joZhb0D/+8CbZY7tXQufngxlW1wJU0REpLHUzRTv1vwfrH7JlJNSYPgkyP1Z4PvdfwXTTofib6F0PXxxCZzyKSR53YlXREQkRGqZiWd718PC2wKvj321fiID0LwLjPwMmnUyrwtnwvJHIxejiIhIEymZiWeLfgdVJaZ82DWQN7rh85p1gOP/BZbHvP76T1C8IjIxioiINJGSmXi1/UvY8G9TTmsPRz1+8PPbDoM+Y025phIWjQ1vfCIiIg5RMhOvlt4dKA+4H1KyDv2eAeMgfd+aM5s+hMLPwxKaiIiIk5TMxKMdC2DrZ6bcoocZ5BuM5HST+Pgtvcf52ERERBymZCYeffdEoNzvjtBmJnUbA5m9TXnbLNg2x9nYREREHKZkJt6UboT1/zTl1DaQd0Vo709Khr53Bl5/95hzsYmIiISBkpl4s+Y1sKtMucf1kNws9GvkjYa0Dqb842TYs8ax8ERERJymZCae2DasfT3wukeQY2V+ypMKvX6775o1sHJC02MTEREJEyUz8WT7HChZacrtR0Dzro2/Vo/rAmNt1v7d7O8kIiIShZTMxJO6rTLdxjTtWmltoNN5plxeCBunNO16IiIiYaJkJl7UVMH6fYvkedIh9+dNv2b3XwbKa/6v6dcTEREJAyUz8WLbLKjcacqdzgJvi6Zfs8MoSO9syps+1I7aIiISlZTMxIsNkwPlzj878HmhSPJAt6tM2a6GDe86c10REREHKZmJB7YNP75nykleyDnTuWt3uSRQ9q9fIyIiEkWUzMSDnQuhdIMptx8Z3D5MwWo5ILAicOFMKNvs3LVFREQcoGQmHmz8b6Cc61AXk59lQZeL972wYf07zl5fRESkiZTMxIPNUwPlnLOcv35tMgOsf9v564uIiDSBkplYV7kLdn5pylmHQ3on5++RdThk9jHlbV9A+Tbn7yEiItJISmZi3ZbPzJYDAB3yw3MPy4LO+xbQw4ZNH4XnPiIiIo2gZCbWbanTxdQxTMkMQM7ZgfKmD8J3HxERkRApmYlltg2bPzHlpBRod2L47tVmGKRkm/Kmj6G6Mnz3EhERCYGSmVi2Zw3s/cGU2x4Pyenhu1eSJ7B+TVWJWXFYREQkCiiZiWV1E4r2I8J/v051upo2qqtJRESig5KZWFZYJ5lpe0L479fxNLCSTVnJjIiIRAklM7HM3zKT5IXWQ8N/v5SW0PY4U96zCvasC/89RUREDsHVZGbmzJmcc8455OTkYFkW7733Xu33fD4fd955JwMGDKB58+bk5ORw1VVXsWnTJvcCjiZlW6BkpSlnHw3JzSJz3w6jAuUtn0bmniIiIgfhajKzd+9eBg4cyDPPPLPf90pLS1m0aBH33XcfixYt4t133+X777/n3HPPdSHSKLTt80C5XQS6mPw6nBooK5kREZEokOzmzc844wzOOOOMBr+XlZVFQUFBvWNPP/00Q4cOZf369XTp0qXB91VUVFBRUVH7uri4GDAtPT6fz6HIqb1m3T8jKWnLDDz7ylXZw7AjFUPGQJK9WVi+Iuwt/6OqsgKsQ+fEbtZVrFFdBU91FTzVVfBUV6EJV32Fcj3Ltm3b0bs3kmVZTJ48mfPPP/+A53z66afk5+eze/duMjMzGzxn3LhxjB8/fr/jEydOJD09jFOXI+yksrG0rFmDjcWH6W9QZbWI2L2Hlj9Mx+p5AExLe5xiz2ERu7eIiCSG0tJSRo8eTVFR0QF/5vvFTDJTXl7O8ccfT58+fXjzzTcPeJ2GWmZyc3PZvn37ISsjVD6fj4KCAkaNGoXX63X02gdVtZfkya2xqMHO6k9V/qLI3RtIWjUBz+JbAKg+4mFqev/ukO9xra5ikOoqeKqr4Kmugqe6Ck246qu4uJg2bdoElcy42s0ULJ/Px6WXXkpNTQ3PPffcQc9NTU0lNTV1v+NerzdsD2U4r92gXV8DZj8mq82xkf/P1uk0WGyKnsJpePrfFfRbI15XMUx1FTzVVfBUV8FTXYXG6foK5VpRPzXb5/Nx8cUXs3btWgoKChxvXYlJO74MlCMxJfunMnpBemdT3jYTqisOfr6IiEgYRXUy409kVq5cyaeffkrr1q3dDik6uJ3MWBa0P8WUq8th54LIxyAiIrKPq91Me/bsYdWqVbWv165dy5IlS8jOziYnJ4cLL7yQRYsW8cEHH1BdXc2WLVsAyM7OJiUlxa2w3bdzvvnT0wyyDncnhnYnwtq/m3LhTGg73J04REQk4bnaMrNgwQIGDRrEoEGDABg7diyDBg3ij3/8Iz/++CPvv/8+P/74I0ceeSQdO3as/Zo9e7abYburfLvZYBIg+yhIcikfbVtnh+7Cme7EICIigsstMyeffDIHm0wVJROtoou/VQbMyr9uyegBaR2gfAts+wJqqtxLrEREJKFF9ZgZacCOOsmMG+Nl/CzLdDUBVJXA7qXuxSIiIglNyUyscXvwb13t1NUkIiLuUzITa/wzh1JaQQuXV95VMiMiIlFAgxxiSdkWKN9qyq2OMl09bso63CRVlbtg2yywa4Lapykh+Yphzeuw8QMo+hp8JZDaGrKHQOfzoctF4EngGXoiIk2gnzyxZFedcSmtBroXh5+VBG337dhdsQOKvnU3nmhk18CKp+C9LrDwZtgyFco2mXFGe9fBhn/DnCvg/cNg3STQoHcRkZApmYklu5cEyq2OdCuK+up2NW2b5V4c0ahyN3yWDwtvAV9R4HhqW8jsDckZgWNlG2H2ZfD5xVC1J+KhiojEMiUzsaRuy0zLKGiZgfqL5W2f414c0aZsMxQcD1v/FzjWbQycuQwu2ApnfwcX7oKR0yHnzMA5G/5N8mcnkVqzM9IRi4jELCUzscQ//TkpBTL7uBuLX6tBJh5QMuNXWQTTToeib8zr1LZwyv9g2GvQckBgrFOSB9qfBCdPgRPeAa/Zd8wq+orh5X8MjI8SEZGDUjITK6rKoPg7U87qFz2DRT2pkD3YlEtWmhWKE1lNNcz6OexeZl437wr5c6DDKQd/X+4FkD/PnA9k2D+SPPNM8KnLSUTkUJTMxIqib8xgUoie8TJ+bYYFyjvmuhdHNPj6gUDXUmobGDEVMroH996sPjByGnZ6F8C00DDnysC/u4iINEjJTKzYtSRQjpbxMn51k5lE7mraOh2+/pMpWx44YTJk9grtGi26UXXCf/GRbl7/+B58/WcnoxQRiTtKZmJF3e0CorllJlGTmaoymPcrYN/U6iP+BO2Ob9y1MvsyP/V2bP9/z6/Hw7YErVcRkSAomYkV0bbGTF3pnSA915R3fGk2nUw03/wZ9qw25bbHQ787m3S5bcmDqOl3j3lhV8Ps0WbhPRER2Y+SmVhg27D7K1NO72xW3Y02/taZqr1mhdtEUvw9fPsXU07ywtAXHFkJuabvPdDmOPNi7zpYck+TrykiEo+UzMSC8i3g223KWYe7GsoBJXJX09J7ocZnyn3vMLPNnJCUDMe9CZ5942dWPgfbE3yAtYhIA5TMxAL/eiUQG8lMIo3v2D7PbEkAkNYeDr/b2eu36AZHPLDvhQ1fXhtInEREBFAyExuKlgfKTv3W77RWgyAp1ZQTpWXGtmHJXYHXA+6H5ObO36f3LaZ+wXQ3rpzg/D1ERGKYkplYEAstM56UwOJ5e1YlxuJ52z6Hwumm3KIHdP9VeO6TlGzG4fh9NQ4qtN2BiIifkplYULdlJrOve3EcSutjAuWd892LI1KWPxoo97/PDP4Nl9ZHQ96Vply5C75+MHz3EhGJMUpmop1tB1pm0jtDSpa78RxM66MD5R1xnszsWgabpphyehfIuyz89xz4Z/A0M+WVz0DxyvDfU0QkBiiZiXblW81v4gCZUTpexq/10EB5x5fuxREJ/qnYAH1/F95WGb/mudD396Zc44Nl94X/niIiMUDJTLSLhcG/fi0Og5RsU94537QqxaPSjfDDJFNObQ3dfxm5e/e9A9LamfL6twPrD4mIJDAlM9EuFgb/+lkWZA8x5fJCKN3gbjzhsvplsyovQI/fhGcG04F4W0C/OjOolt0fuXuLiEQpJTPRLpZaZiD+x83U+GDVi6ZsJUGPayMfQ4/roVlHU/5xMuxcFPkYRESiSLLbAcgh1GuZibFkZud86PJz92IJhx/fh7JNptzpXDOOJdKSm0G/e2DhTeb1sj/CyR9EPo5YUeMz0+i3fQG7vzYraleVmha1tPbQsr9Z9LHtCWaJARGJOUpmol3xvpaZZjmQ0tLVUIKSXbdlJg4HAa98LlDueYN7cfT4tRmEXLrBzKrasQBaD3Evnmi0Zw189yT88BZUHGTdo/Vvmz+TM6DLhebfVXUpElPUzRTNKnaYL4DMPu7GEqz0HJN4AexcCHaNu/E4qWQVbP3MlFv0gA4j3YvFkwqH19l4cvkj7sUSbcoLYe4v4b+94PunD57I1FVVAmtehU+Ohv+NjM9uUpE4pWQmmpXUWUcks7d7cYTKP0XbV2x2lI4Xa/8eKPf4tSM7YzfJYb+AtA6mvOFdKPrO1XBcZ9uw5u/wQV9Y83+BQdqeNOhyCRz7Opz1LVy8By6rgYv3wlnfwLA3Ie8K8GYGrrX1M/hkKMy+MjFWsxaJcUpmolndRCCjl3txhOqn42bigV0Da98wZSvJ/PBzmycN+ty274Vdf+2bRFNdBnOvhrljoHLfVg/eLBjwAPxsMxw/CQ67CrL6mLEylgXJ6WYcWrfL4bg34GebzLYRLboHrrvuTZjSD9b/252/l4gERclMNCuJ0WQmOw5nNG37HPauM+UOo0x3WjToeT14W5ry2jdgb5xOhz+ItJqdJH92Eqx9PXCw62Vwzvcw4L7gx5olNzez085aDkOegZRW5njFNvj8Iph/A1RXOB6/iDSdkploVjeZyYyhZKbu4Ml4SWbW1PlB2e0q9+L4KW8m9LrRlO0q+O4xd+OJtD2rOb78bqzdS8zr5OYwfBIMnxhYXDBUnhRTp2d9A53PDxxfOQEKjoe9PzQ1ahFxmJKZaFa8wvxpJUPzPFdDCUlKKzNAFmDXYqiudDeepqoqhfX/MuXkjPo/4KJB71sCezateilxxngUfUfytBE0t7ea183zIH8edL3Emes36wgnvAvHvGK69AB2LjBjaeIlSReJE0pmopVdExgAnNEdkmJsFr1/3ExNBRR97W4sTfXj+2amC0CXi8xYi2iS1ha6/8qUq0vNDJ54t2cdfHYqVvkWAOzMfjDqC2jp8CrZlgXdr4H8OYGxNOWF8OnJ5rkQkaigZCZalW40gxohtsbL+MXTSsDr/xkod4uCgb8N6ft704IH8P0zULXX3XjCqWwLfDYKyjYCsDvpMKpGfBbecUytjoTTvoR2J5rX1aUw62eB1aBFxFUx9ut+AonVwb9+2T+Z0ZR3jXuxNIVvD2z+yJTT2kHbE92N50CadzGDXte9YWbzrH4Fet/sdlTOqyqDmefBnlUA2Bm9mFN9L6f6NzgNp9RsGDEV5v7CbDRq18CX15luyD63hv/+saZ8m1k4c8d808pc+gOUbjK/pFWXm1mB3ixIyYLmXU3XdGZvs7RDVn9I8rj9N5AYomQmWsXq4F+/7EHmw8quMavTxqqNH5gPXoDcn0f3B2y/200yA/DtY9DzN5DkdTcmJ9k2fPnrwMrS6blUnfgRldMjuHO4JxWO+wekd4Zv/585tug20xLW/97IxRGNbNuMKfrxP7Dx/eB2dK/cCXuBXUvqH09uYbaYyDkLOp8DLQ4LR8QSR5TMRKtYXWPGL7m52eV791dmzExVqdsRNU7dLqYuF7kXRzBaDoCcM2HTh1C6Hn74p1lDJV58+1dY9w9TTm4OJ/0X0nOBCCYzYJL0I/9iBoN/tW/X8mV/MF1PRzxoxtkkkvLtZkHJ1S9D8bcHPze1tUlUklLNooa+YqjcZWbi1VW1B7YUmK9Ft5rPkrzRkHelO/uhSdRTMhOtYr2bCSB7iElm7GqsomVuRxM6X0lsdDHV1e9Ok8yAWUQvb3R8/HDd8j9Yclfg9bC/Q6uB4PO5E49lwYA/msHgi283x755yHSDHfVYfNT5oexdbxLM1S8HWi9rWWbcXNvjofUxJtFu3rXhwfM11Wb8U8lK2LUUts+B7bMDG7qC2XB36b2w9A/Q/hQ47Gqzj5YnNax/RYkdSmailb9lJrm5mSIai1ofbfa6AaydC4AYayqOpS4mv7YnmB8eO+bB7mWw+WPIOcPtqJqmbCvMvgKwzesB4yD3AjcjCuj7e/Ckw4J9a/2seAJqKmHIU+5vdxEupZtMi9Sa1/ZvUWl3omk96XQ2NOsQ3PWSPGbMV/Mugf3ObNskMBv/a762z9l3sg1b/2e+Fv8OelxnFo6M1c/IcKquhPKtZtHH8m1QucN8ntX4zJddDcnNzPObnG7WrErrYL5Ss2Pu+VUyE42qK2HvWlPO6BW7v+VlBxbPs3YtJOaSGf/aMhD9XUx+lmVaZ2bt+2G//NHYTmbsGphzFeybgk3H06D/fe7G9FO9bjDr0Mz7FWDDymfB9sHRE2LuB8JBVe0144SW/8V0qfklN4fu15qkwqnxfZYFLfubr8PvNlPx175hVnnes9qcU74Vvn7AtIh1ucgMeG9zrDP3jyVVpaZFa+cCKFpuBseXrDQtZ/5fAEJlJZsEMaMHZPQ0g7MzeprtP1p0j8pf7FxNZmbOnMlf//pXFi5cyObNm5k8eTLnn39+7fdt22b8+PG8+OKL7Nq1i2OOOYZnn32Www93eC2JaLN3bWCTvFjtYgJoeYQZgFrj29cyEyMJAZgPiFjrYvLrfJ55bkq+h8IZsH0etDnG7aga59u/wpapppzWwXQvRWOC0P0aSEoxe0PZNWbKdk0lDH05Kj/4Q7ZhMiz4bf2uH2+WSSB632LGwoRTizyzNUX/P8C2WWb5gQ3vms9Juwp+eMt8tT4Get8KXX4eX4Pf6/IVm27XTR+ZFqvi5eaZc5JdBaUbzNfWafW/50k3n+2tjjRdva2OhOZ9nb1/I7iazOzdu5eBAwdy9dVX8/Of/3y/7//lL3/h8ccf57XXXqNXr148+OCDjBo1ihUrVpCRkeFCxBFSsjpQzujhXhxN5Uk1D/3OhVDyPcnpMTQIeMungS6mTufF1g8kKwn63m5m/oAZO3PCO+7G1Bg7FphxEgBYcNybjd+iIBK6XWF+o51zhfkhu+Y108o67PXYW/TSr3w7LLzJTEX3s5Kh5w1mzFC4k5ifsizTldXuRLMP2arnTeJYsW/V6x3zYPZlsLgT9Pqt2d0+0jE6zd/ltulDk8Bs+3z/7r2f8maZlpT0XLOoZmpbSG1jupMsr0n0rCTzGVdValraKneZ1q6yLaYltHQ9VOzY/9rVpbBjrvnaJ5kkRlidoepkc28XuPo/7IwzzuCMMxpuArdtmyeffJJ7772XCy4wTeavv/467du3Z+LEiVx33XWRDDWy9tRJZuru4BuLso+GnQuxsMmqWX3o86PFxjqru3Y+1704GqvblfDVH6Fss/mtuniFWcMjVlSX72vl2NdCefg9gfEU0SzvUvOD4otL97UYTDR/Hvdm7LUUrH8H5v/GjLnwyzkTjno8Op6l5rkw8M+m23HdW7DiSTNODMyA4qV3m26oblea1qOsfq6GGxJfyb7Wlw9NC3Hpjw2fZ3nMmjzZg82eeC0HmlbZ1NbODE+o3AXFK023Vcn3ZkLH7qWwZ039MKghmTJXV0eP2l8X1q5dy5YtW8jPz689lpqaykknncTs2bMPmMxUVFRQURHY2ba4uBgAn8+Hz+GZD/7rOX3dpOJV+NsBqpp1wXZrxoYDrJaDah+yljWrHa+rsLBrSP7xv1iA7WlGVesTIzprxpnnKomkHjfh+eoewKbmm79SPWSCI/FFQtKyP+ApWg6A3XIQVX3uafDfIFz/B5uk47lYx/0Tz5xLsWoqYf0/qamuoPrYf5iuKJcEXVcV2/AsuoWkH/9de8j2tqJ60BPYXS4zPySjqb7xQJcrIPdyrG0zSVr5FNamD7CwzQJ9q16EVS9S034UNT1vwu6Qf8iuyog/V7YNxd+QtOUTrM2fYG3/Astu+N5288Oo6Xg6dofTsNue1HACUXWIlptgWS0ga5D5qstXjFX0FdbupVi7l8LOhezck0l2mH7GBhWqbduNHCHkLMuy6o2ZmT17NsOHD2fjxo3k5ASWKb/22mv54Ycf+OSTTxq8zrhx4xg/fvx+xydOnEh6epTtqXMAQ8sfomO1WRjsk2YvUZ7U1uWIGi+zZh0jym4FYKNnOAvSbnc3oCC0ql7BieV3ArDZM5Qv0+5xOaLGSbb3kl/6a7yUUk0yBc1epCIpAivlNlGr6u84ofweLGqoJpkZzR6jJKmr22GFrF3VQoZWPIIH84G8xTOE+al3UGO5l9AcSseq2QyseIFUimqPbfYMZWnKb6hIauViZKFJr9nMYb4pdKn6H17K6n2vxOrEOu9pbPSc4OrfKdkuo031MtpXL6Rd9SLS7YY3iK3Gy3ZPfwo9R7HVcxR7rZzonBRi247HVVpayujRoykqKiIzM/Og50Z9MrNp0yY6dgxMu/v1r3/Nhg0b+Pjjjxu8TkMtM7m5uWzfvv2QlREqn89HQUEBo0aNwut1rgk5+ZNBWMXfYCelUHVBcXQOeAxWTRXJ77XGqi5jr9Ue67w1jtZVOCR99Qc83/0FgKohL2F3GxPR+zv5XCUtuxvPiscAqO79e2qOeMiJEMOnqpTkgqOx9phNVqsHPEhNnzsOeHq4/g86xdr6PzxfXIC1b5+1mvb5VA//V2CX8wg6aF2VF+JZfAtJPwbGVtkp2VQPehI795Lo/OEZDF8xSeteJ2nls1h763eN2CRhtx9JTZdLsXPOhpRAYhOW58quxtq1BGtrAdbWT7G2zzl460uH07A7nobd9uTo29z2J8L1/7C4uJg2bdoElcxEbTdThw5mjYItW7bUS2YKCwtp3779Ad+XmppKaur+Cyl5vd6wfdg5em3brp2WbbXohjcl1heF8kKrQbB9Ns3trfhqivF6g1x/wi2bp+wrWCR3ORdc+iHpyHPVdyysfBpqKvGseRHPgD+YvXCi1bLxsC+RofVQPIffiSeIwbPh/P/dJJ1Ph5M/ghlnQdVekrZOJemLn8GJ/wFvC1dCqldXtm2WIFhwY2AQLUDnn2Ed/RzJwa4VE628raHfWOhzC2z6AL57EgqnA2ach7W1gKStBeYXxjbDoOMZZlG+jP7m7U15rqrKzOSHHXPNrKPCGQ0PqAWzInK7k8wyCjlnYGX0whODCaTT/w9DuVbUJjPdunWjQ4cOFBQUMGiQ6a+rrKxkxowZPProoy5HF0blWwNrODSPsXVZDiR7iFnRE7B2LYIWZ7oc0EGUrDYzB8B8uEXz7JlgpOeYAZCrXzFTOle9aPZwikbbvjCDOMGs23JsDM8Cqqv9STDiE5h2BlSVwNbP4NMT9m3H0Nm9uMq2mOnWG+rMdEttDUOehS4Xx25rTEOSPGbJgs7nQdF3ZluMdW/C3nXm+3aNef62fQFAspXMSVZXPPPfg6ye5rM4PcfM1PFm1S45QY0PfLvN53b5VjMwtvg7KPrWDJo92KyjFt3Nukk5Z0D7EWa9Hmk0Vz8p9uzZw6pVq2pfr127liVLlpCdnU2XLl249dZbeeihh+jZsyc9e/bkoYceIj09ndGjR7sYdZjVHSUeL5urtQ7soG3tWgC5UZzMbPxvoByLs5ga0vd2WP1/gG1WqO19c/QtA19dDvN+Se0iX0c8CFl9XA3JUW2HwykFMO008BWZjRU/GWpaaOr8/4gI24Y1r5sNMit3BY7nXgBDnoNmB275jgtZfWDgn+CIB8wvWevfMTOGir+rPcWyq2hpr4Z1Ds7A9GaZGXkd8qHjqPj5fI8SriYzCxYsYMSIEbWvx44dC8CYMWN47bXXuOOOOygrK+OGG26oXTRv6tSp8b3GTN1kJiPGp2X71V0JeOdCFwMJQt0p2Z3iJJnJ7A2dz4cfJ5up2uvehO6/dDuq+r56wEwfh8DCZ/GmzTGQPxtmnGP+n5dthk9PNC0hh10dkZaQZjWFeGadA1unBg6mttnXGnNRfLXGHIplmSSz7XDgcdizFjZPhR1zsbfPg+LvzIyoxkhKNashtx5qWnhbHwuZfWJrvaoY42oyc/LJJ3Ow8ceWZTFu3DjGjRsXuaDcVm+NmTjJ3DN7YSdnYFWV7NvWIEpVFkHhLFNu0d18+MSLfneaZAbMqrqHXR09A8t3LjIL+4Fpvj/mlfj90M/qB/nzzHYT22YFWqS2FMDRz4dvPFNVGUnL/8IpZQ+TVBaYIEHX0TD4SbOwWqJr0Q16Xgc9r6PK52PqlHfIP+4wvOXrzedy+TbTVesrMusfWcnmefVm7NvTqD2kd4LMvtA8L36f4SgVBx3ScSYeu5msJOxWR2Ftm4FVttH8RhqNG8Nt/SzQx51zZnz9ltrmGLNqauFM0wLy4/uQe77bUZkxB/N+WWdxvPugZZxvV5LWxnQ5LbzZjGECs8Luti9gyNNmXIdT7BozJmbxHXj840MAmnUye0d1Pse5e8WZKquZmbzgHep2KBKEKPnVTGrVTWaad3MvDofZ2YMDL3YscC+Qg9lcZ7p/x9PdiyNc+t4ZKC9/1IydcNu3/8+MHwFoOcC0ICUCTyoMfQGO/2dg+ffSDTDzfJh+dqBOGqumGn74J3w4ED6/uHagaw1JVPe4Ec76RomMxBUlM9HGn8yktXNt6mY42K3qJDM7ozCZsW3YtC+ZSUo1M1DiTc4ZZulzMNNFt81yN56i7+CrfQtcWkmme8kTvQvKhUWXi+CMJWZQqN+mKfDRIJhxHmz8EGpCWM117w9m/NH73eCLS6Do69pv1bQbyfRmT1Iz6Inonp4v0gjqZoomVWWBXWljfU+mn7DrDAJmx3z3AjmQ4u/MxmpgumPicZqkZUG/O2DOVeb18kfN39UNdo3pXqrZN36jz+8iP6snWrTIgxEfm/VeFo01+wqBGYy+8X0zHqPDKPNvldnbTOf2pJlducu3mWd35yLY+qnZO+enWh8LA+6nus0plHz0UUT/aiKRomQmmuxbLA+In/Eyful5VJBBKiWmZSYMS183Sbx3Mfl1vRSW/sEkbps+hG2zoe1xkY/j++dq1x6iRQ8YMC7yMUQTy4KuF5vlAFa9BMsfCfxiU74F1r1hvoK+XpJZAK73LdDh1CjcT0nEWepmiibxOPjXz7LY7elhyhXbAq0g0WJTnWQmJ46TmSQv9L838HrJXZEfO7NnHSy9K/D6mJejfrn2iPGkQe+b4Lx1cOJ7Zkq9J9i6scy09iMehHPXwskfmPVMoumXBpEwUctMNCmJw2nZdexO6kH76sXmxY4F0DxKNg+sKjVLjQOk55qplfHssGvg28eg5HszbmbTFOh0dmTuXVMNc8dA1V7zusf18Tk+qamSvIEVa6srzRinXUvMqrLlhaZ7zkqG1LbQvIuZdZM9WFOsJWEpmYkm9Vpm4mvMDJhkptbO+dDl5+4FU1fhjMDYjY6nx/9vsknJMPAh+PxC83rJ3aZLIhLrYqx40kwPB5PMDorjrUmc4kkx42XcGt8kEgPUzRRN6q4D0TzPrSjCpl4yE03TsxOli6mu3AvM6qRgZrys+b/w33P317D0nn0vLLP3ktfZnexFJDEpmYkme38wf1rJ0CzH3VjCoDypNXbavsXydi4wM1qigX/wr+WB9iPdjSVSLAuO/Evg9dK7D7yjrxOqK80sqppK87rPWHUviYhjlMxEE38yk945bpfCtlsdZQq+ovpjhNyyZ40ZOwLQ5rjEWn+j/UlmdhOYRGbJ3eG715I7Yde+8VJZ/WDgg+G7l4gkHCUz0cJXbLaSh+gZGBsG9dabiYbF8zZ/EignShdTXYMeg+R9G7eufhm2zXH+HuvfMWNlAJJSYNibZtaOiIhDlMxEC3+rDMR3MlN3JeBoWDxvU4KsL3Mg6TlwxL5VeLFNV5B/ppETSlbBvGsCrwc/CdmDnLu+iAhKZqJHoiQz2VG0rUF1JWz9nymntYNWR7oajmt63WTWJwHYswoW/c6Z61buhpnnmVZHMF1aPa535toiInUomYkWCZLMmHUx9v39di0y6464ZfsXgVaIDqeZVVMTUVIyDHsjsDjbqhdg/b+bds0aH8y6EIqWm9eZfWDoi/E/7V1EXJGgn95RqF4yk+daGBGRvW8Pnqq9Zl8ZtyTilOwDyewJRz0eeD1nTON3bq6phjm/CLR6pbaBkz4Ab0ZToxQRaZCSmWiRKC0zAK3rDgJ2cdxM7X5MltnIL9H1uBbyrjDl6lKYfnboM878K/z+MNG8Tko1y/JnxN8ikCISPZTMRIu6yUx6rntxREJ2nd2R3Vo8r3QT7F62L54hWgYeTBfQMS8Fxs+UbYT/jYDiFcG9v7IIZpwN6/5hXid54fh/Qtvh4YlXRGQfJTPRwp/MNOsInlR3Ywm37KMCZbdaZhJ9SvaBeNLgpPch63DzunQDfDLUTK8+mMJZ8PGQQGtXkheO/7fZBVpEJMyUzESD6nIo32LK6XHexQSQ0hIyepnyrqVmVlGkbU7wKdkHk9YORn4GLY8wr33FZh+naafDpo+gqswcr66ELZ/CrJ/DpyeamVAAqa1hRIESGRGJGG00GQ32bgiU4328jF/2ELPybk2F2RuobmtNuNVUwZYCU/a2DOxRJAFp7WDUFzDvl7D+n+bY5k/Ml+UxCWnlrv23pGgzDI57My53fReR6KWWmWhQmkCDf/1a1xk3E+n1ZnbMNz+IATqOMlOTZX/eFjB8kvmqO47LrjbbH9RNZJp1hCHPwKmzlMiISMTpUzwaJNJMJr+62xrsmG9m0kSKupiCZ1nQ9RKzy/bGD2Dj+6Zr0FcEKa3M2Jqcs0yXkrYoEBGXKJmJBgmZzAwyi9TZNZFvmamXzORH9t6xKskLuT8zXyIiUUbdTNEgEZOZ5OaQ2c+Ud38VGFQabhU7AntCZfU3O5SLiEhMUzITDRIxmYHAuBm7GnYvjcw9NxcAtilrSraISFxQMhMN/MlMSqvEWvK93riZLyNzT42XERGJO0pm3FZTDaU/mnIitcoAtDk2UN72RfjvZ9uBxfI86dD2+PDfU0REwk7JjNvKt4JdZcrxvo3BT7U8ApJbmPK2z02yEU67lwUWJ2w/Iv5XWhYRSRBKZtzmb5WBxEtmkpIDrTNlm+qPHQoHdTGJiMQlJTNuK62z+m8izqyp29Wz7fPw3mtTnWRGg39FROKGkhm3JXLLDPwkmQnjuBlfSSBZatEdMnqE714iIhJRSmbcVlY3mUnAlpnWx5i9fiC8LTNbPwuMTVIXk4hIXGlUMuPz+diwYQMrVqxg586dTseUWEoTPJnxtoBWg0y56OvAnklO2/hBoKwuJhGRuBJ0MrNnzx5eeOEFTj75ZLKyssjLy6Nfv360bduWrl278utf/5r58+eHM9b4VHfMTLNO7sXhprbDA+Vts52/vl0TSGY8adD+FOfvISIirgkqmXniiSfIy8vjpZde4pRTTuHdd99lyZIlrFixgjlz5nD//fdTVVXFqFGjOP3001m5cmW4444f/paZ1DaQ3MzdWNwS7nEzOxfVmZJ9KiSnO38PERFxTVAbTc6ePZtp06YxYMCABr8/dOhQrrnmGp5//nleeeUVZsyYQc+ePR0NNC7ZNVC60ZQTsYvJr24yUzjd+etv/G+g3Pkc568vIiKuCiqZ+de//hXUxVJTU7nhhhuaFFBCKS8MDEptlsDJTLMOkNkHir8z2xr4isGb6dz16yYzOWc5d10REYkKms3kpkRfY6au9iPNn3Y1FM507rqlP8Kuxabc6ihIT9BxSSIicSyolpm6ysvLefrpp5k2bRqFhYXU1NTU+/6iRYscCy7u1Z3J1DwB15ipq8OpsPJZU97yP+h0tjPXrTuLqZO6mERE4lHIycw111xDQUEBF154IUOHDsWyrHDElRjqJjOJ3M0E0P5ksJLMOKItnzp3XY2XERGJeyEnM1OmTOHDDz9k+PDhhz65iaqqqhg3bhz/+Mc/2LJlCx07duQXv/gFf/jDH0hKioMeskRfY6aulJbQajDsnG/WmynbCs3aN+2alUWBxKhZjulmEhGRuBNyMtOpUycyMjLCEct+Hn30UZ5//nlef/11Dj/8cBYsWMDVV19NVlYWt9xyS0RiCCuNmamvw6kmmQGzYm/eZU273sb/Qk2lKef+HNSKKCISl0JOZh577DHuvPNOnn/+ebp27RqOmGrNmTOH8847j7POMjNQ8vLyeOutt1iwYMEB31NRUUFFRUXt6+LiYsCsWuzz+RyNz3+9xl7Xs3dD7Qhsn7c9OBxfNAmmrqw2J5HMwwDUbPqE6k4XNumenh/erq3fqpyfYcdI/Tb1uUokqqvgqa6Cp7oKTbjqK5TrWbZt26FcfNu2bVx88cXMnDmT9PR0vF5vve87ub3BI488wvPPP8/UqVPp1asXS5cuJT8/nyeffJLLLmv4t/Zx48Yxfvz4/Y5PnDiR9PToWizt1NLraG5vpZIMPmr+htvhuC7JruSM0itJpoJysvgk/f8C+zaFKNku5fTSMXjwUW614pNmLzf6WiIiEnmlpaWMHj2aoqIiMjMPvlxHyMnMqaeeyvr16/nlL39J+/bt9xsAPGbMmNAjPgDbtrnnnnt49NFH8Xg8VFdX8+c//5m77777gO9pqGUmNzeX7du3H7IyQuXz+SgoKGDUqFH7JXWHZNeQ/G4mVk0ldtYAqvIXOhpbtAm2rjxfXEjSpvcBqDplJnbrYxt1P2v9WyTPM89idfffUHPU3xp1HTc06blKMKqr4Kmugqe6Ck246qu4uJg2bdoElcyE3M00e/Zs5syZw8CBAxsdYLDefvtt3nzzTSZOnMjhhx/OkiVLuPXWW8nJyTlg0pSamkpqaup+x71eb9geykZdu2xr7XgOq3luwvyHOWRd5Z4P+5KZ5C0fQocTGnejje/WFj15l+CJwfoN5zMbb1RXwVNdBU91FRqn6yuUa4WczPTp04eysrJQ39Yot99+O3fddReXXnopAAMGDOCHH37g4YcfdrQFyBVlmsnUoE5nARZgw8b34ciHQ79G+TbYOMWU0zrU3y5BRETiTsjzmx955BF+97vfMX36dHbs2EFxcXG9LyeVlpbuNwXb4/Hst1BfTKo3LTvBF8yrK60dtBlmykXLoWRV6Nf44a3ANhHdroQkjZUREYlnIbfMnH766QCMHDmy3nHbtrEsi+rqamciA8455xz+/Oc/06VLFw4//HAWL17M448/zjXXXOPYPVyjNWYOrPN5sH22Ka//Jxx+T2jvX/N6oNwtxlvwRETkkEJOZqZNmxaOOBr09NNPc99993HDDTdQWFhITk4O1113HX/84x8jFkPYaI2ZA+tyMSy505TXvgH97g5+jZhdy2DXvi01so+GloeHJ0YREYkaISczJ510UjjiaFBGRgZPPvkkTz75ZMTuGTGlmwLlZtr8sJ4WedD2BNg2y+ykvWsxZAe5eu/KCYHyYWqVERFJBCGPmXn11Vf517/+td/xf/3rX7z++usNvEMaVFY3mclxL45o1e2KQHnt34N7T+WuwLnJLSDvioOfLyIicaFRA4DbtGmz3/F27drx0EMPORJUQvAnM5508Dq7/k1c6HIRJO2bYr/mNfDtOfR7Vr0M1aWmfNgvICUrXNGJiEgUCTmZ+eGHH+jWrdt+x7t27cr69esdCSoh+JOZZjnaM6ghKa0g73JT9hXBukOskFxVCt89Hnjd66bwxSYiIlEl5GSmXbt2LFu2bL/jS5cupXXr1o4EFfeqSs0PaIB0dTEdUO+bA+VvH4Oag+zTsXIClG8x5dwLILNXeGMTEZGoEXIyc+mll3LzzTczbdo0qqurqa6u5rPPPuOWW26pXdxODqFsc6Cc1tG9OKJdq4HQft8SAHtWw+qXGz6vfBt84+/itGDA/ntziYhI/Ao5mXnwwQc55phjGDlyJM2aNaNZs2bk5+dzyimnaMxMsDT4N3gD6zxTX42D8u37n7NoLFTu2+A073Jo2T8ioYmISHQIeWp2SkoKb7/9Ng8++CBLliyhWbNmDBgwgK5du4YjvvhUN5lRN9PBtRlqBgOv/xeUF8Lcq+HE9wKr+q56Eda9acreljDor25FKiIiLgk5mfHr2bMnPXv2dDKWxFG3m0ktM4c2+G+wdRpUbIdNH8DnF0G/u2DLVFhWZwHFIU9Bsw7uxSkiIq4IqpvpkUceobS0NKgLzps3jylTpjQpqLhXr5tJY2YOqVlHGPYmWPty7x8nw9RjYNl9gG2O9b7V7MMkIiIJJ6hkZvny5XTp0oXf/OY3fPTRR2zbtq32e1VVVSxbtoznnnuO4447jksvvZTMTK2bclAaMxO6nNPgxMlmynZdVhL0vw+Oerzh94mISNwLqpvp73//O8uWLePZZ5/l8ssvp6ioCI/HQ2pqam2LzaBBg7j22msZM2YMqampYQ065imZaZxOZ8PZK8wiesXLzUywvNEa8CsikuCCHjNzxBFH8MILL/D888+zbNky1q1bR1lZGW3atOHII49scFVgOQB/MpPcArwZ7sYSa9LaQr/b3Y5CRESiSMgDgC3LYuDAgQwcODAc8SQG/wBgjZcRERFpspDXmZEm8u0BX7Epq4tJRESkyZTMRJqmZYuIiDhKyUykafCviIiIo5TMRFq9lhmNmREREWkqJTORppYZERERR4U8m2nv3r088sgj/O9//6OwsJCampp631+zZo1jwcUlJTMiIiKOCjmZ+dWvfsWMGTO48sor6dixI5ZlhSOu+KVkRkRExFEhJzMfffQRU6ZMYfjw4eGIJ/5pzIyIiIijQh4z06pVK7Kzs8MRS2KoXf03A7wt3I1FREQkDoSczPzpT3/ij3/8Y9C7aMtP+JOZdHUxiYiIOCHkbqbHHnuM1atX0759e/Ly8vB6vfW+v2jRIseCizu+EqjaY8oaLyMiIuKIkJOZ888/PwxhJIi6g3/TNF5GRETECSEnM/fff3844kgMdQf/qptJRETEESEnM34LFy7k22+/xbIs+vXrx6BBg5yMKz7VTWbSOrgXh4iISBwJOZkpLCzk0ksvZfr06bRs2RLbtikqKmLEiBFMmjSJtm3bhiPO+FC+NVDWtGwRERFHhDyb6aabbqK4uJhvvvmGnTt3smvXLr7++muKi4u5+eabwxFj/CjfEiintXcvDhERkTgScsvMxx9/zKeffkrfvn1rj/Xr149nn32W/Px8R4OLO3VbZtTNJCIi4oiQW2Zqamr2m44N4PV699unSX6irE7LTDMlMyIiIk4IOZk55ZRTuOWWW9i0KTDNeOPGjdx2222MHDnS0eDijr+byUqGlFbuxiIiIhInQk5mnnnmGUpKSsjLy6N79+706NGDbt26UVJSwtNPPx2OGOOHv5sprT1YIVe9iIiINCDkMTO5ubksWrSIgoICvvvuO2zbpl+/fpx66qnhiC9+1FRDeaEpq4tJRETEMY1eZ2bUqFGMGjXKyVjiW+UOsKtNWTOZREREHBNUMvPUU09x7bXXkpaWxlNPPXXQczU9+wA0k0lERCQsgkpmnnjiCS6//HLS0tJ44oknDnieZVlKZg5EM5lERETCIqhkZu3atQ2WJQRaME9ERCQsQp5S88ADD1BaWrrf8bKyMh544AFHgopL6mYSEREJi5CTmfHjx7Nnz579jpeWljJ+/HhHgopL6mYSEREJi5CTGdu2sSxrv+NLly4lOzvbkaDq2rhxI1dccQWtW7cmPT2dI488koULFzp+n7BTN5OIiEhYBD01u1WrVliWhWVZ9OrVq15CU11dzZ49e7j++usdDW7Xrl0MHz6cESNG8NFHH9GuXTtWr15Ny5YtHb1PRKibSUREJCyCTmaefPJJbNvmmmuuYfz48WRlZdV+LyUlhby8PIYNG+ZocI8++ii5ubm8+uqrtcfy8vIcvUfE+LuZklLBm+luLCIiInEk6GRmzJgxAHTr1o3jjjuuwc0mnfb+++9z2mmncdFFFzFjxgw6derEDTfcwK9//esDvqeiooKKiora18XFxQD4fD58Pp+j8fmvF8x1k8u3YgF2WgeqqqocjSMWhFJXiU51FTzVVfBUV8FTXYUmXPUVyvUs27btQ51UXFxMZmZmbflg/Oc5IS0tDYCxY8dy0UUX8eWXX3LrrbfywgsvcNVVVzX4nnHjxjU4EHnixImkp6c7FlsoLLuKc0ovwsJmZ1IvZjX7iytxiIiIxIrS0lJGjx5NUVHRIXOLoJIZj8fD5s2badeuHUlJSQ0OAPYPDK6urm585D+RkpLCkCFDmD17du2xm2++mfnz5zNnzpwG39NQy0xubi7bt293NNECkzUWFBQwatSog7dUlW3C+0EeADU5Z1M9/F1H44gFQdeVqK5CoLoKnuoqeKqr0ISrvoqLi2nTpk1QyUxQ3UyfffZZ7UyladOmNT3CIHXs2JF+/frVO9a3b1/eeeedA74nNTWV1NTU/Y57vd6wPZSHvHbJztpiUnoOSQn8nyOc/w7xRnUVPNVV8FRXwVNdhcbp+grlWkElMyeddFKD5XAbPnw4K1asqHfs+++/p2vXrhGLwRH1pmVrJpOIiIiTQl5n5uOPP+bzzz+vff3ss89y5JFHMnr0aHbt2uVocLfddhtz587loYceYtWqVUycOJEXX3yRG2+80dH7hF2Z1pgREREJl5CTmdtvv712EPBXX33F2LFjOfPMM1mzZg1jx451NLijjz6ayZMn89Zbb9G/f3/+9Kc/8eSTT3L55Zc7ep+wq7vGjFb/FRERcVTQU7P91q5dWzuO5Z133uGcc87hoYceYtGiRZx55pmOB3j22Wdz9tlnO37diFI3k4iISNiE3DKTkpJSu9Hkp59+Sn5+PgDZ2dmHnLadsNTNJCIiEjYht8wcf/zxjB07luHDh/Pll1/y9ttvA2ZgbufOnR0PMC7U28pAyYyIiIiTQm6ZeeaZZ0hOTubf//43EyZMoFOnTgB89NFHnH766Y4HGBf83UzJLcDbwt1YRERE4kzILTNdunThgw8+2O/4E0884UhAccnfMqNWGREREceFnMyA2SX7vffe49tvv8WyLPr27ct5552Hx+NxOr7YV10JlfumrCuZERERcVzIycyqVas488wz2bhxI71798a2bb7//ntyc3OZMmUK3bt3D0ecsatiW6CsZEZERMRxIY+Zufnmm+nevTsbNmxg0aJFLF68mPXr19OtWzduvvnmcMQY28oLA+W0tu7FISIiEqdCbpmZMWMGc+fOrd2rCaB169Y88sgjDB8+3NHg4kLdZCa1nXtxiIiIxKmQW2ZSU1MpKSnZ7/iePXtISUlxJKi4Uq+bScmMiIiI00JOZs4++2yuvfZa5s2bh23b2LbN3Llzuf766zn33HPDEWNsq9fNpGRGRETEaSEnM0899RTdu3dn2LBhpKWlkZaWxvDhw+nRowd/+9vfwhFjbKuo282kMTMiIiJOC3nMTMuWLfnPf/7DypUr+fbbbwHo168fPXr0cDy4uKCWGRERkbBq1DozAD179qxNYCzLciyguKNkRkREJKxC7mYCeOWVV+jfv39tN1P//v15+eWXnY4tPtQOALYgpbWroYiIiMSjkFtm7rvvPp544gluuukmhg0bBsCcOXO47bbbWLduHQ8++KDjQcY0f8tMahtI0grJIiIiTgs5mZkwYQIvvfQSl112We2xc889lyOOOIKbbrpJycxP+ZMZLZgnIiISFiF3M1VXVzNkyJD9jg8ePJiqqipHgoobVXuhutSUtWCeiIhIWISczFxxxRVMmDBhv+Mvvvgil19+uSNBxQ0N/hUREQm7Rs1meuWVV5g6dSrHHnssAHPnzmXDhg1cddVVjB07tva8xx9/3JkoY1W5Vv8VEREJt5CTma+//pqjjjoKgNWrVwPQtm1b2rZty9dff117nqZrowXzREREIiDkZGbatGnhiCM+qZtJREQk7Bq1zowEScmMiIhI2CmZCae6yYxmM4mIiISFkplwqtAAYBERkXBTMhNO9bqZNABYREQkHJTMhJN/NpOVDN6WroYiIiISr5TMhFPtVgbtQFPVRUREwkLJTLjYdmDMjMbLiIiIhI2SmXDxFUGNz5Q1k0lERCRslMyEiwb/ioiIRISSmXDRGjMiIiIRoWQmXCq0+q+IiEgkKJkJF21lICIiEhFKZsKlvM7qv9oxW0REJGyUzISLuplEREQiQslMuKibSUREJCKUzISLkhkREZGIUDITLv5uJk8zSG7ubiwiIiJxTMlMuFRsN3+mtnE3DhERkTinZCYcbBsqdpqykhkREZGwiqlk5uGHH8ayLG699Va3Qzk4XzHYVaac2trdWEREROJczCQz8+fP58UXX+SII45wO5RDq9wRKKcomREREQmnmEhm9uzZw+WXX85LL71Eq1at3A7n0CrqJDNqmREREQmrZLcDCMaNN97IWWedxamnnsqDDz540HMrKiqoqKiofV1cXAyAz+fD5/M5Gpf/ej+9rlW6tbZiq5NbUuPwfWPRgepK9qe6Cp7qKniqq+CprkITrvoK5XpRn8xMmjSJRYsWMX/+/KDOf/jhhxk/fvx+x6dOnUp6errT4QFQUFBQ73XnqhkM3ldevnora9Z/GJb7xqKf1pUcmOoqeKqr4Kmugqe6Co3T9VVaWhr0uVGdzGzYsIFbbrmFqVOnkpaWFtR77r77bsaOHVv7uri4mNzcXPLz88nMzHQ0Pp/PR0FBAaNGjcLr9dYeT1q5GpaYct8jT6BP1zMdvW8sOlBdyf5UV8FTXQVPdRU81VVowlVf/p6VYER1MrNw4UIKCwsZPHhw7bHq6mpmzpzJM888Q0VFBR6Pp957UlNTSU1N3e9aXq83bA/lfteu2l1bTE5vD/rPUCuc/w7xRnUVPNVV8FRXwVNdhcbp+grlWlGdzIwcOZKvvvqq3rGrr76aPn36cOedd+6XyEQNDQAWERGJmKhOZjIyMujfv3+9Y82bN6d169b7HY8qSmZEREQiJiamZsccrTMjIiISMVHdMtOQ6dOnux3CoflbZqxk8Do76FhERETqU8tMONRuMpkNluVuLCIiInFOyUw4+FtmtMmkiIhI2CmZcVp1OVTvW+hH42VERETCTsmM0zSTSUREJKKUzDhNyYyIiEhEKZlxmqZli4iIRJSSGaepZUZERCSilMw4zT8tG5TMiIiIRICSGafVa5nR1GwREZFwUzLjtAqNmREREYkkJTNOq9SYGRERkUhSMuM0DQAWERGJKCUzTqvXzZTtXhwiIiIJQsmM0/zdTN4sSIq5TclFRERijpIZp9XumK2ZTCIiIpGgZMZJNdVQuduUNV5GREQkIpTMOKlyF2CbsqZli4iIRISSGSdpWraIiEjEKZlxkqZli4iIRJySGSdp9V8REZGIUzLjJG0yKSIiEnFKZpxUqU0mRUREIk3JjJM0ZkZERCTilMw4ScmMiIhIxCmZcVKlBgCLiIhEmpIZJ6llRkREJOKUzDjJn8x40iA53d1YREREEoSSGSf5p2ari0lERCRilMw4xbYDY2Y0LVtERCRilMw4pWoP1PhMWeNlREREIkbJjFM0+FdERMQVSmacomnZIiIirlAy4xS1zIiIiLhCyYxTtMmkiIiIK5TMOKVCm0yKiIi4QcmMUyo0ZkZERMQNSmacUqkxMyIiIm5QMuMUDQAWERFxhZIZpyiZERERcYWSGaf4u5msJPBmuRuLiIhIAlEy45TaTSazTUIjIiIiERHVP3Uffvhhjj76aDIyMmjXrh3nn38+K1ascDushlVok0kRERE3RHUyM2PGDG688Ubmzp1LQUEBVVVV5Ofns3fvXrdDq6+m0mw0CRovIyIiEmHJbgdwMB9//HG916+++irt2rVj4cKFnHjiiS5F1QCtMSMiIuKaqE5mfqqoqAiA7OzsA55TUVFBRUVF7evi4mIAfD4fPp/P0Xj816sq3Yp337EabyuqHb5PPPDXldP/BvFIdRU81VXwVFfBU12FJlz1Fcr1LNu2bUfvHia2bXPeeeexa9cuZs2adcDzxo0bx/jx4/c7PnHiRNLT08MSW+vqrzm+/A8ArPSez/KUX4TlPiIiIomitLSU0aNHU1RURGZm5kHPjZlk5sYbb2TKlCl8/vnndO7c+YDnNdQyk5uby/bt2w9ZGaHy+XwUFBRw+uGlpH45GoDqAQ9S0+cOR+8TD/x1NWrUKLxe76HfkMBUV8FTXQVPdRU81VVowlVfxcXFtGnTJqhkJia6mW666Sbef/99Zs6cedBEBiA1NZXU1NT9jnu93rA9lJ7q3YFys/Z49PAfUDj/HeKN6ip4qqvgqa6Cp7oKjdP1Fcq1ojqZsW2bm266icmTJzN9+nS6devmdkgNsrT6r4iIiGuiOpm58cYbmThxIv/5z3/IyMhgy5YtAGRlZdGsWTOXo6ujcmegrGRGREQkoqJ6nZkJEyZQVFTEySefTMeOHWu/3n77bbdDq8eq1NRsERERt0R1y0yMjE3WJpMiIiIuiuqWmZhRt5sp5cBr4IiIiIjzlMw4wPJvMpmcAZ4Ud4MRERFJMEpmnOBvmdEmkyIiIhGnZKap7Bqo3GXKGi8jIiIScUpmmsjLXixqzAslMyIiIhGnZKaJUuySOi+UzIiIiESakpkmqpfMqGVGREQk4pTMNJGSGREREXcpmWmiFLu4zgslMyIiIpGmZKaJUqjbMqOp2SIiIpGmZKaJ1M0kIiLiLiUzTeRVMiMiIuIqJTNNpJYZERERdymZaSKtMyMiIuIuJTNNVDsAOCkFkpu7G4yIiEgCUjLTRLVTs1PbgGW5G4yIiEgCUjLTFLYd6GbSeBkRERFXKJlpiupSPPhMWeNlREREXKFkpikqdwTKapkRERFxhZKZpqhQMiMiIuI2JTNNYFXuDLxQN5OIiIgrlMw0RcX2QFktMyIiIq5QMtME9VpmtMmkiIiIK5TMNIUGAIuIiLhOyUxTVGjMjIiIiNuUzDSBpZYZERER1ymZaYp6Y2aUzIiIiLhByUxT7FtnxsYCb0t3YxEREUlQSmaaoLabKaUVJHncDUZERCRBKZlpCv8KwBr8KyIi4holM41V48OqKgbA1ngZERER1yiZaax607Kz3YtDREQkwSmZaay607LVzSQiIuIaJTONVWdfJnUziYiIuEfJTGOVbw2UU9u5F4eIiEiCUzLTWGWBZMZOUzIjIiLiFiUzjVVRGCintXcvDhERkQSnZKax6nQz2epmEhERcY2SmcYq/i5QbpbjXhwiIiIJTslMY1SVwfZ5AOyxOqqbSURExEUxkcw899xzdOvWjbS0NAYPHsysWbPcDSi5GZy/gapj/8GKlEvdjUVERCTBRX0y8/bbb3Prrbdy7733snjxYk444QTOOOMM1q9f725gaW2xcy/ix+ST3I1DREQkwUV9MvP444/zy1/+kl/96lf07duXJ598ktzcXCZMmOB2aCIiIhIFkt0O4GAqKytZuHAhd911V73j+fn5zJ49u8H3VFRUUFFRUfu6uNhsBunz+fD5fI7G57+e09eNR6qr4Kmugqe6Cp7qKniqq9CEq75CuZ5l27bt6N0dtGnTJjp16sQXX3zBcccdV3v8oYce4vXXX2fFihX7vWfcuHGMHz9+v+MTJ04kPT09rPGKiIiIM0pLSxk9ejRFRUVkZmYe9Nyobpnxsyyr3mvbtvc75nf33XczduzY2tfFxcXk5uaSn59/yMoIlc/no6CggFGjRuH1eh29drxRXQVPdRU81VXwVFfBU12FJlz15e9ZCUZUJzNt2rTB4/GwZcuWescLCwtp377h6dCpqamkpqbud9zr9YbtoQznteON6ip4qqvgqa6Cp7oKnuoqNE7XVyjXiuoBwCkpKQwePJiCgoJ6xwsKCup1O4mIiEjiiuqWGYCxY8dy5ZVXMmTIEIYNG8aLL77I+vXruf76690OTURERKJA1Cczl1xyCTt27OCBBx5g8+bN9O/fnw8//JCuXbu6HZqIiIhEgahPZgBuuOEGbrjhBrfDEBERkSgU1WNmRERERA5FyYyIiIjENCUzIiIiEtOUzIiIiEhMi4kBwE3h360hlJUEg+Xz+SgtLaW4uFgLKx2C6ip4qqvgqa6Cp7oKnuoqNOGqL//P7WB2XYr7ZKakpASA3NxclyMRERGRUJWUlJCVlXXQc6J6o0kn1NTUsGnTJjIyMg64n1Nj+fd92rBhg+P7PsUb1VXwVFfBU10FT3UVPNVVaMJVX7ZtU1JSQk5ODklJBx8VE/ctM0lJSXTu3Dms98jMzNQDHyTVVfBUV8FTXQVPdRU81VVowlFfh2qR8dMAYBEREYlpSmZEREQkpimZaYLU1FTuv/9+UlNT3Q4l6qmugqe6Cp7qKniqq+CprkITDfUV9wOARUREJL6pZUZERERimpIZERERiWlKZkRERCSmKZkRERGRmKZk5iCee+45unXrRlpaGoMHD2bWrFkHPX/GjBkMHjyYtLQ0DjvsMJ5//vkIRRodQqmv6dOnY1nWfl/fffddBCOOvJkzZ3LOOeeQk5ODZVm89957h3xPIj9XodZXoj5XDz/8MEcffTQZGRm0a9eO888/nxUrVhzyfYn4bDWmrhL1uQKYMGECRxxxRO2CeMOGDeOjjz466HvceK6UzBzA22+/za233sq9997L4sWLOeGEEzjjjDNYv359g+evXbuWM888kxNOOIHFixdzzz33cPPNN/POO+9EOHJ3hFpffitWrGDz5s21Xz179oxQxO7Yu3cvAwcO5Jlnngnq/ER/rkKtL79Ee65mzJjBjTfeyNy5cykoKKCqqor8/Hz27t17wPck6rPVmLryS7TnCqBz58488sgjLFiwgAULFnDKKadw3nnn8c033zR4vmvPlS0NGjp0qH399dfXO9anTx/7rrvuavD8O+64w+7Tp0+9Y9ddd5197LHHhi3GaBJqfU2bNs0G7F27dkUguugE2JMnTz7oOYn+XNUVTH3puTIKCwttwJ4xY8YBz9GzZQRTV3qu6mvVqpX98ssvN/g9t54rtcw0oLKykoULF5Kfn1/veH5+PrNnz27wPXPmzNnv/NNOO40FCxbg8/nCFms0aEx9+Q0aNIiOHTsycuRIpk2bFs4wY1IiP1dNkejPVVFREQDZ2dkHPEfPlhFMXfkl+nNVXV3NpEmT2Lt3L8OGDWvwHLeeKyUzDdi+fTvV1dW0b9++3vH27duzZcuWBt+zZcuWBs+vqqpi+/btYYs1GjSmvjp27MiLL77IO++8w7vvvkvv3r0ZOXIkM2fOjETIMSORn6vG0HNldhoeO3Ysxx9/PP379z/geXq2gq+rRH+uvvrqK1q0aEFqairXX389kydPpl+/fg2e69ZzFfe7ZjeFZVn1Xtu2vd+xQ53f0PF4FUp99e7dm969e9e+HjZsGBs2bOD//b//x4knnhjWOGNNoj9XodBzBb/97W9ZtmwZn3/++SHPTfRnK9i6SvTnqnfv3ixZsoTdu3fzzjvvMGbMGGbMmHHAhMaN50otMw1o06YNHo9nv1aFwsLC/TJOvw4dOjR4fnJyMq1btw5brNGgMfXVkGOPPZaVK1c6HV5MS+TnyimJ9FzddNNNvP/++0ybNo3OnTsf9NxEf7ZCqauGJNJzlZKSQo8ePRgyZAgPP/wwAwcO5G9/+1uD57r1XCmZaUBKSgqDBw+moKCg3vGCggKOO+64Bt8zbNiw/c6fOnUqQ4YMwev1hi3WaNCY+mrI4sWL6dixo9PhxbREfq6ckgjPlW3b/Pa3v+Xdd9/ls88+o1u3bod8T6I+W42pq4YkwnN1ILZtU1FR0eD3XHuuwjq8OIZNmjTJ9nq99iuvvGIvX77cvvXWW+3mzZvb69ats23btu+66y77yiuvrD1/zZo1dnp6un3bbbfZy5cvt1955RXb6/Xa//73v936K0RUqPX1xBNP2JMnT7a///57++uvv7bvuusuG7Dfeecdt/4KEVFSUmIvXrzYXrx4sQ3Yjz/+uL148WL7hx9+sG1bz9VPhVpfifpc/eY3v7GzsrLs6dOn25s3b679Ki0trT1Hz5bRmLpK1OfKtm377rvvtmfOnGmvXbvWXrZsmX3PPffYSUlJ9tSpU23bjp7nSsnMQTz77LN2165d7ZSUFPuoo46qN3VvzJgx9kknnVTv/OnTp9uDBg2yU1JS7Ly8PHvChAkRjthdodTXo48+anfv3t1OS0uzW7VqZR9//PH2lClTXIg6svxTPH/6NWbMGNu29Vz9VKj1lajPVUN1BNivvvpq7Tl6tozG1FWiPle2bdvXXHNN7ed627Zt7ZEjR9YmMrYdPc+VZdv7RuaIiIiIxCCNmREREZGYpmRGREREYpqSGREREYlpSmZEREQkpimZERERkZimZEZERERimpIZERERiWlKZkRERCSmKZkRkYiZPn06lmWxe/duV+7/2Wef0adPH2pqag557gcffMCgQYOCOldE3KVkRkTC4uSTT+bWW2+td+y4445j8+bNZGVluRLTHXfcwb333ktS0qE/+s4++2wsy2LixIkRiExEmkLJjIhETEpKCh06dMCyrIjfe/bs2axcuZKLLroo6PdcffXVPP3002GMSkScoGRGRBz3i1/8ghkzZvC3v/0Ny7KwLIt169bt18302muv0bJlSz744AN69+5Neno6F154IXv37uX1118nLy+PVq1acdNNN1FdXV17/crKSu644w46depE8+bNOeaYY5g+ffpBY5o0aRL5+fmkpaXVHlu6dCkjRowgIyODzMxMBg8ezIIFC2q/f+655/Lll1+yZs0aR+tHRJyV7HYAIhJ//va3v/H999/Tv39/HnjgAQDatm3LunXr9ju3tLSUp556ikmTJlFSUsIFF1zABRdcQMuWLfnwww9Zs2YNP//5zzn++OO55JJLANNism7dOiZNmkROTg6TJ0/m9NNP56uvvqJnz54NxjRz5kwuu+yyescuv/xyBg0axIQJE/B4PCxZsgSv11v7/a5du9KuXTtmzZrFYYcd5lDtiIjTlMyIiOOysrJISUkhPT2dDh06HPRcn8/HhAkT6N69OwAXXnghb7zxBlu3bqVFixb069ePESNGMG3aNC655BJWr17NW2+9xY8//khOTg4Av//97/n444959dVXeeihhxq8z7p162rP91u/fj233347ffr0AWgwEerUqVODSZiIRA8lMyLiqvT09NpEBqB9+/bk5eXRokWLescKCwsBWLRoEbZt06tXr3rXqaiooHXr1ge8T1lZWb0uJoCxY8fyq1/9ijfeeINTTz2Viy66qF4sAM2aNaO0tLTRfz8RCT8lMyLiqrrdOgCWZTV4zD9FuqamBo/Hw8KFC/F4PPXOq5sA/VSbNm3YtWtXvWPjxo1j9OjRTJkyhY8++oj777+fSZMm8bOf/az2nJ07d9K2bdtG/d1EJDKUzIhIWKSkpNQbtOuUQYMGUV1dTWFhISeccEJI71u+fPl+x3v16kWvXr247bbbuOyyy3j11Vdrk5ny8nJWr17NoEGDHItfRJyn2UwiEhZ5eXnMmzePdevWsX37dscWn+vVqxeXX345V111Fe+++y5r165l/vz5PProo3z44YcHfN9pp53G559/Xvu6rKyM3/72t0yfPp0ffviBL774gvnz59O3b9/ac+bOnUtqairDhg1zJHYRCQ8lMyISFr///e/xeDz069ePtm3bsn79eseu/eqrr3LVVVfxu9/9jt69e3Puuecyb948cnNzD/ieK664guXLl7NixQoAPB4PO3bs4KqrrqJXr15cfPHFnHHGGYwfP772PW+99RaXX3456enpjsUuIs6zbNu23Q5CRCQS7rjjDoqKinjhhRcOee62bdvo06cPCxYsoFu3bhGITkQaSy0zIpIw7r33Xrp27RrUWJ61a9fy3HPPKZERiQFqmREREZGYppYZERERiWlKZkRERCSmKZkRERGRmKZkRkRERGKakhkRERGJaUpmREREJKYpmREREZGYpmRGREREYpqSGREREYlp/x+OZsKixUzt1AAAAABJRU5ErkJggg==\n",
"text/plain": [
"